A parameter-free metadynamics 
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We present a unique derivation of metadynamics. The starting point for the derivation is an on- 
the-fiy reweighting scheme but through an approximation we recover the standard metadynamics and 
the well-tempered metadynamics in a general form while never appealing to the extended Lagrangian 
framework. This work leads to a more robust understanding of the error in the computed free energy 
than what has been obtained previously. Moreover, a formula for the exact free energy is introduced. 
The formula can be used to post-process any existing well-tempered metadynamics data allowing 
one, in principle, to obtain an exact free energy regardless the metadynamics parameters. 



The last decade has seen the introduction of a number 
of adaptive biasing techniques developed for free energy 
computation^ - — In practice, these techniques accumulate 
a biasing potential or a biasing force during trajectory 
evolution. The goal is to specify a biasing force or poten- 
tial that will effectively flatten the free energy landscape 
and enhance sampling. Here we focus on a particular 
adaptive biasing potential (ABP) method called metady- 
namics, introduced in references [H and §■ Metadynamics 
is an ABP method that has been widely used for chemi- 
cal, solid state and biological systems. 7 

One area of metadynamics that can still be improved 
upon is that of the dependence of the error in the com- 
puted free energy on the metadynamics parameters^— 
Metadynamics requires one to specify several system de- 
pendent parameters. The parameters are the energy rate, 
which is the product of the Gaussian height and the de- 
position period, and the Gaussian widths The error as- 
sociated with these is mostly well understood, with one 
exception: The Gaussian width. At present, the an- 
alytic error estimates suggest that, for fixed computa- 
tional time, the error decreases as the Gaussian width 
increases^ while numerical experiments reveal that the 
error will actually increase in the limit of large Gaussian 
widths and that the empirically derived error estimate 
will break down£. 

Below, we derive metadynamics in a novel way. Origi- 
nally, the goal was to propose an adaptive biasing scheme 
with on-the-fly reweighting (as in references and 
then, by an approximation, show that the reweighting 
factors could be removed. The problem with the on- 
the-fly reweighting methods is that they are difficult to 
study formally because of the time-dependent reweight- 
ing. The hope was to introduce an ABP method that 
was free from these reweighting factors, allowing it to be 
studied rigorously like metadynamics and adaptive bi- 
asing force methods^ However, this process leads one 
directly to metadynamics, both the well-tempered and 
standard forms. This derivation accentuates an inter- 
pretation that leads to an understanding of the error 
incurred by choosing a finite Gaussian width and ulti- 
mately affords an exact expression for the free energy, 
independent of the Gaussian width. In other words, a for- 
mula follows from this interpretation that is exact even 



for finite Gaussian width. This formula could be used 
to post-process any existing well-tempered metadynam- 
ics data. We also stress that the biased dynamics can 
be cast in a form consistent with the dynamics studied 
in reference Hit The biasing force can be expressed as an 
average taken over replicas of the system rather than an 
average over time. It is expected that because of this, it 
should be possible to obtain convergence results for the 
well-tempered metadynamics without assuming "instan- 
taneous equilibration" of the dyanamics. 

Below we present the derivation by first introducing 
the on-the-fly reweighting strategy and then introducing 
our approximations. Once metadynamics is recovered we 
present an error estimate and an exact formula for free 
energy. 

Let x be a single configuration in the n-dimensional 
configuration space X of some interesting dynamical sys- 
tem. For this system, assume N collective variables 
(CVs) have been specified and further suppose that the 
CVs are good descriptors of the interesting features of A". 
The CV space is fi. Let £ denote a point in f2. Following 
reference y, let us define the mollified free energy (up to 
an arbitrary constant £) 



where 



Z- 1 f x 6 a (ax)-0e- mx) dx, 



or 



(1) 



(2) 



Let Z absorb the normalization of S a . In the limit 
a — > this is the exact free energy. When a is finite, 
the exact free energy can in principal be recovered via 
deconvolution^ and we will later make use of this. 

In practice, equation ([TJ would be computed via trajec- 
tory averages. Let xt be a trajectory solving the following 
Langevin equation, 



-fx t dt - VV{x t )dt + y/2fp- l dB t 



where we assume it to be ergodic with respect to the 
Boltzmann density on A", / is the Langevin friction co- 
efficient and j3 is the inverse temperature. The random 
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force dB t is given by the increments of a Brownian mo- 
tion. 

Denning the population at £ 

g(t,t)= [ 6 a {t{x,)-{)da (3) 



and 



Z t = f g(£,t)dt, 
Jn 



(4) 



the free energy may be computed from a long trajectory 
xt as 

Ce -^(M = Zr 1 g(U), (5) 

where at t = 



(6) 



Now, we propose the following on-the-fly reweighting 
scheme to compute equation (|SJ) via an adaptive biasing 
potential. In this case we propose the biasing potential 



(7) 



where b and c are scalars. b controls the strength of the 
bias and c > is a coupling parameter with units of 
inverse time. With the bias defined this way, the initial 
conditions for the biased dynamics reduce to those of the 
unbiased case and < V&. 

Using this bias potential, one can write the mollified 
free energy as time averages over Langevin trajectories 
driven by V + V& with on-the-fly reweighting, 



where 



and 



Z t = f g(t,t)dt. 
Jn 



(8) 



(9) 



(10) 



It is natural for the development of the ABP method to 
terminate here; these equations could be implemented 
as-is. See reference for example. 

In this case one can go a bit further, however, by in- 
troducing an approximation to equation ©. The idea 
behind the following manipulations is to derive an ex- 
pression for g(£,,t) in which the factor exp[/3Vf,] does 
not appear. Once obtained, we recover mctadynamics. 
Using a first-order expansion Vb(£,(x t ),t) w Vb(£,t) + 
Vfr(£,t)(£(xt) — £) and an expansion of e x , we obtain 



Jo 

= / Sa(a*s) - 0^ VbiLs) ds + A(a). (11) 



The last term is roughly the error between the popula- 
tions computed by equation © and the population com- 
puted with the first term in the last line of equation , 

A(a) = 4^ t d e 5 a (ax s )-0n&s)e^' s) ds. (12) 



Ignoring A (a), the final line in equation (fTT|) can be 
seen as an integral solution to the differential equation 



dt 



S a (Z(xt))-0e mm , 



which can be solved via separation of variables 

g(£, t) = - c ^ [c{l - b) J* 8 a (t(x a ) - Ods + I)' " - l^j 

(13) 

To avoid the possibility of complex valued g, we restrict 
c > and b < 1. We can now express g without evaluat- 
ing exp[/3Vh]- 

Equations (J7]) and (TT31 combine to give the bias po- 
tential 

V6fet) = /3- 1 r ^ln[c(l-6) f 5 a (t(x.)-t)da + l]. 

(14) 

In the limit b — > 1 the bias reduces to 

V b (Z,t) = /T x c I 6 a (£(x a ) - £)ds. (15) 



These bias potentials are the well-tempered^ and stan- 
dard metadynamics^ bias potentials, respectively. In- 
deed, making the substitutions uj = j3~ 1 cb and AT = 
/3 _1 &/(1 — b) in equation ([T4| . we find the well-tempered 
metadynamics, 

V b (U) = AT In (JL J - Ods + lj . (16) 



Equation (|15p is the standard metadynamics with an en- 
ergy rate of /3 _1 c. We have derived metadynamics via 
the approximations in equation (1111) having started with 
an on-the-fly reweighting scheme. 

To derive an exact expression for the free energy when 
< & < 1, it will be useful to work with the biasing force 
and with averages taken over independent replicas of the 
dynamics, rather than averages in time. The gradient of 
equation (fl"4"|) is 



(17) 



l + c(l-&) / 5 a (Z(x s )-0ds 
Jo 



where 



( AA(i(.o) .,, = 2 (6(*«) & 5 a Mxt)-t). (18) 



a* 
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In the long-time limit, 



(1-6) / 5 a (£(x a )-t)da 
Jo 



(19) 



Notice that in the long-time limit c vanishes. This im- 
plies that for well-tempered metadynamics the Gaussian 
height will not impact the long-time accuracy of the com- 
putation. 

Noting the i—th replica of the biased dynamics with 
xt(i), the replica density can be defined 

1 M 



»=i 



Each replica is an independent solution to the biased 
Langevin-equation 



dx t = -fx t dt - V{V{x t ) + V b ^{x t )))dt + ^/2f(i- 1 dB t . 

(20) 

Equation (fT9| can be cast as the following replica average 



(21) 



^ (l-b) [ 6 a (t(x) - M(x,t)dx 

J x 

Assuming that when t — > oo, the method converges 
and ip(x, oo) oc exp[— j3(V(x) + Vb(£(x), oo))} so that 

r 1 / d^ix) - Oe- p{V(x)+Vbii(xhoc)) dx 

J X 



x 



r 1 / dt5 a ($-t)e 
Jn 



S a (Z-Z)e 



-0(A{£)+V b {t,oo)) 



dn a (£,oo) 



(22) 



where we have defined the mollified free energy in the bi- 
ased ensemble fi a . For finite a, the biasing force is thus a 
rescaling of the mollified mean force in the biased ensem- 
ble and V b — —bfi a /(l — b) up to an additive constant. 
This interpretation of the biasing force in well-tempered 
metadynamics leads directly to an exact expression for 
the free energy. 

Consider the exact mean force in the biased ensemble, 



= -(V e A(0 + %«H(f,oo)) 

= -v fi .A(r) + r ^.|i a (r,oo). 



(23) 



The exact free energy to an unimportant constant is 
therefore, 



A(£) =^,oo) + — -/i a (£,oo). (24) 

What makes this expression useful is that /i can be ob- 
tained by a deconvolution of \i a and the known Gaussian 

If one ignores the deconvolution and assumes /i w [i a , 



A(0-~lv b (C,oo) 



(25) 



where the error associated with this approximation is 
roughly, 

/T 1 In (l + ^ [Wfa oo)) 2 ~ oo)] ) • (26) 

We have used Vf, = — b /i a /(l — 6) in equation (|2"5"|) . This 
error estimate was also given in reference H] (see equation 
(21) there). This estimate is found by making a Taylor 
series expansion of exp[— /3/i(£, oo)] in 



and keeping terms up to the second moment of 5 a . No- 
tice that the error in well-tempered metadynamics is re- 
lated to a convolution of the configurational density in 
the biased ensemble. In reference |j the error was due to 
a convolution of the configurational density in the unbi- 
ased ensemble. 

In practice one needs the histograms 



(28) 



and the exact free energy can be computed with 

MZ, *) = -P' 1 In {5~ x * t) x h b (li, t)) (29) 

where we use S^ 1 * h to indicate a deconvolution. 

By now it has been demonstrated that metadynamics 
is a powerful computational toot, so we use a very sim- 
ple model here to demonstrate that the formula given in 
equation (|29|) can be used to compute an accurate free en- 
ergy even for large a. The Richardson-Lucy schem o 12 i 13 
was used for the deconvolution, as described in reference 
4. 

We apply the method with £(x) — x and A(x) — 



V(x) = 



+ 0.25. The dynamics are given by 



equation O where / = l/2di, k B T = .25/10 (one 
tenth of the barrier height), the particle mass is unity. 
The Langevin integrator from reference [l4] was used to 
evolve the trajectory. The coordinate x is discretized 
from x = —2 to x = 2 into 400 bins. At t = the initial 
phase point is (x, x) — (\/0.5, 0). 
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Here we implement the above stated grid-based meta- 
dynamics for a single trajectory with b = 0.8 and c = 
1/dt. The histogram h(£,t) defined above and its deriva- 
tive 

h'fot)= [ d £i 8 a (Z(*s)-0ds (30) 
Jo 

are computed on the grid of values £, where at each 
timestep the trajectory makes a contribution to all grid 
points. The gradient of the bias potential is given by 
equation (|17j) . At the end of the simulation we apply 
equation (|29[) to remove the impact of finite a. 

We evolve for 10 7 dynamical steps and compute the 
error e = — A(£i,i)|d£/4 where i runs over all 

bins such that < /3 _1 and 4 is the length of the 

grid. We adopt this condition from reference |8j. The 
error is shown in figure [1] where A is either the estimate 
of A with or without deconvolution. The two are clearly 
labeled in figure [1] and we find that the deconvolution 
makes a significant improvement to the computed free 
energy. The grid spacing should satisfy d£ « a so that 
the 6 a are well represented on the grid. If the grid is too 



coarse, one can expect an increase in error as a decreases. 

In conclusion, we have presented a derivation of meta- 
dynamics that leads to an understanding of the error as- 
sociated with finite a and a formula for removing this 
error. We have demonstrated that the formula in equa- 
tion can indeed be used to correct the computed free 
energy even when a is large. Equation (|29[) can be used 
to post porcess any existing well-tempered metadynam- 
ics data to remove the blurring related to using a finite 
a. In hindsight, it appears that both equation ([55)1 and 
the error in equation (l26l) should follow straight from the 
presentation of metadynamics in reference [l|. One only 
needs to notice that the biasing force is related to the 
mollified mean force in the biased ensemble. 
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FIG. 1: Error as a function of a with and without deconvo- 
lution. 



